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LOCATING  AN  APPROPRIATE  SADDLEPOINT 
FOR  H-DIHENSIONAL  PROBABILITY  INTEGRALS 

INTRODUCTION 

The  M-dimensional  (MD)  moment  generating  function  (MGF)  for  M 
dependent  random  variables  (RVs)  of  interest  can  sometimes  be 
found  in  closed  form;  for  example,  see  references  1  and  2. 
However,  this  joint  MGF  is  usually  not  the  desired  end  product  of 
an  analysis.  Instead,  its  MD  inverse  Laplace  transform,  namely, 
the  joint  probability  density  function  (PDF),  might  be  the 
quantity  of  interest.  Alternatively,  the  joint  MGF  might  be 
weighted  by  an  MD  function  prior  to  performing  the  MD  inverse 
Laplace  transform,  thereby  yielding  the  joint  cumulative 
distribution  function,  or  the  joint  exceedance  distribution 
function,  or  some  other  statistics  of  the  RVs  of  interest;  see 
reference  2,  equation  (18)  and  examples  1  through  7  for  several 
informative  illustrations  of  this  general  approach. 

Since  the  joint  MGF  frequently  exists  only  in  strips  of 
analyticity  in  its  MD  domain,  the  MD  Bromwich  contour  (BC)  for 
returning  to  the  probability  domain  must  initially  be  taken  in 
this  region  of  existence  of  the  joint  MGF.  However,  at  least  for 
one-dimensional  problems,  M  =  1,  the  MGF  can  be  analytically 
continued  outside  this  strip,  and  the  original  BC  can  be  moved  to 
a  more  attractive  region  of  the  complex  plane,  taking  advantage 
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of  possible  saddlepoints  (SPs)  and/or  paths  of  steepest  descent 
of  the  analytic  continuation  of  the  MGF.  in  fact,  for  many 
one-dimensional  examples,  an  exact  result  for  the  first-order  PDF 
can  be  obtained  by  moving  the  BC  to  an  appropriate  valley  at 
infinity,  while  accounting  for  any  residues  of  poles  or  essential 
singularities  traversed  during  the  contour  movement. 

For  M  dimensions,  the  situation  is  much  more  complicated. 

The  analytic  difficulties  associated  with  multiple  complex 
domains  and  the  movement  of  multiple  contours  discourages  and 
often  prevents  any  useful  improvement  in  trying  to  move  the  BCs 
out  of  their  original  MD  domain  of  definition.  Instead,  a  more 
fruitful  approach  is  to  try  to  find  a  single  dominant  SP  of  the 
pertinent  integrand  in  the  original  MD  strip  of  analyticity  and 
develop  a  saddlepoint  approximation  (SPA)  (plus  a  correction  term 
if  possible)  about  this  point  in  MD  space.  However,  it  should  be 
noted  that  there  is  a  question  about  the  existence  of  an  SP  in 
the  original  strip  of  analyticity;  this  problem  will  be 
illustrated  in  this  report  later. 

Even  if  a  single  dominant  SP  exists  in  the  MD  strip  of 
analyticity,  there  remains  the  numerical  problem  of  specifying 
and  quantifying  its  location  in  complex  MD  space,  it  will  be 
shown  that  this  MD  search  procedure  can  be  restricted  to  MD  real 
space  on  a  real  scalar  function  that  has  a  monotonic,  bowl-like 
behavior,  thereby  guaranteeing  a  single  minimum  in  MD  space. 


2 


LOCATING  AN  APPROPRIATE  SADDLEPOINT  IN  MOMENT  GENERATING 
SPACE  FOR  EVALUATION  OF  A  PROBABILITY  DENSITY  FUNCTION 

Let  MD  real  vector  z  =  [z^  Scalar  function  p(z)  is 

an  MD  joint  PDF;  that  is,  p(z)  is  real,  nonnegative,  and  has  unit 
volume  in  MD  z  space.  Let  MD  complex  vector  X  =  +  i  in 

terms  of  its  real  and  imaginary  MD  component  vectors,  The 

MD  Laplace  transform  of  p{z)  is  denoted  as  the  scalar  joint  MGF 

/j(\)  =  J  dz  p(z)  exp(X^  z)  .  (1) 

For  X  purely  real,  that  is,  X^  =  0,  it  is  presumed  that  integral 

(1)  converges  for  real  vector  X^  in  an  MD  real  region  R  that 

includes  the  origin,  X^  =  0,  as  an  interior  point  in  MD  X^  space. 
The  joint  MGF  /y(X^)  is  obviously  positive  real  for  all  X^  c  R^. 

It  follows  that  if  //(X^.)  exists,  then  MGF  fj{\)  =  //(X^.  +  i  X^ ) 

also  exists  for  all  X^,  as  seen  from  the  relation 

|//(X^  +  i  X^)|  =  jj  dz  p(z)  exp^X^  z  +  i  xT  zj  | 

<  J  dz  p(z)  exp^X^  zj  =  //(Xj,)  for  all  X^  .  (2) 

This  MD  region  of  complex  X  is  the  region  of  analyticity  of  joint 
MGF  //(X)  and  is  denoted  by  ROA(//);  that  is,  ROA(//)  is  the  MD  set 
of  complex  X  values  such  that  Re(X)  e  R^.  Joint  MGF  //(X)  does 
not  exist  outside  of  ROA(/t/),  although  its  analytic  continuation 
(AC)  may.  A  one-dimensional  example  is  furnished  by  PDF  p(z)  * 

H  exp(-|z|)  for  all  z;  the  corresponding  MGF  is  //(X)  =  1/(1-X^) 

for  -1  <  Xj.  <  1.  The  AC  of  fj{\)  is  =  1/(1-X^)  for  all 
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^  /  ±1/  at  which 


points  has  poles. 


Equation  (2)  also  indicates  that,  for  a  fixed  real  vector 
e  R^,  the  maximum  value  of  |^(X)|  anywhere  on  the  M  BCs, 

^  ~  ^  ^i'  realized  as  all  the  contours  simultaneously 

cross  the  M  real  axes  in  X  space.  An  alternative  interpretation 
IS  that  the  real  axes  in  X  space  are  mountain  ridges  of  joint  MGF 
U(X)|;  any  perpendicular  departure  from  the  real  axes  in  any  or 

all  of  the  M  complex  planes  {X^}  of  vector  X,  leads  to  a  decrease 
in  |//(X)(  when  X  e  ROA(//). 

The  joint  PDF  p(z)  at  real  MD  field  point  z  can  be  recovered 
from  the  joint  MGF  ^(X)  according  to  inverse  Laplace  transform 

“  Tihr”  I  exp(-x'r  z)  ,  (3) 

where  MD  contour  C  is  a  set  of  vertical  BCs  that  stay  in  the 
ROA(/i).  The  integrand  in  equation  (3)  is  denoted  as 

Y(X,z)  -  //(X)  exp{-X^  z)  . 

The  exponential  can  be  expressed  as  exp(-xj  z)  exp(-i  xT  z), 
which  retains  a  constant  magnitude  as  vector  X^  alone  is  varied. 
Therefore,  integrand  T(X,z)  in  equations  (4)  and  (3)  realizes  its 
maximum  magnitude  on  the  BCs  at  the  real  points  where  the 
individual  contours  all  cross  the  real  axes  in  X  space. 

Therefore,  the  major  contribution  to  integral  (3)  is  expected  to 
occur  in  the  neighborhood  of  the  points  where  the  M  BCs  cross  the 
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real  axes  in  X  space. 


Since  Cauchy's  integral  theorem  states  that  the  value  of 
integral  (3)  is  the  same  for  any  set  of  vertical  BCs  in  the 
ROA(jt/)/  it  is  possible  to  choose,  for  the  crossing  point,  that 
MD  real  point  in  space  where  the  positive  real  integrand 
Y(Xj.,z)  is  a  global  minimum;  this  MD  real  point  is  identified  as 
X  =  X(z).  All  other  axis  crossings  encounter  a  larger  peak 
magnitude  of  integrand  Y(X,z),  and  therefore  must  also  encounter 
more  oscillations  in  the  integrand  away  from  the  crossing  point 
in  order  to  compensate  and  maintain  the  identical  integral  value 
for  PDF  (3)  at  that  z  value.  Selection  of  the  crossing  point  X(z) 
as  the  minimum  of  real  integrand  Y{X^,z)  essentially  minimizes 
the  oscillations  of  complex  integrand  Y(X,z)  on  the  vertical  BCs 
in  the  neighborhood  of  the  crossing  point  and  partially  realizes 
the  steepest  decay  of  the  magnitude  of  the  integrand 
Y(X  +  i  X^,z).  However,  the  selection  of  the  location  of  the 
minimum  of  Y(X^,z)  as  the  crossing  point  X  to  be  used  for  the 
BCs  is  not  mandatory;  it  is  taken  mainly  for  convenience  and 
simplicity.  An  alternative  for  a  one-dimensional  example 
involving  a  branch  point  (BP)  at  the  border  of  is  given  in 
appendix  A. 

This  particular  real  axes  crossing  point,  X(z),  is  the 
location  of  a  maximum  of  |Y(X,z)|  along  every  BC,  and  is  the 
location  of  a  simultaneous  minimum  of  Y(X^,z)  along  every  real 
axis  in  X  space.  If  all  the  slopes  of  Y(X^,z)  are  zero  at 
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X  -  X(2),  this  crossing  point  is  an  SP  in  MD  X  space.  Thus, 

3  ^ 

ax  ^  =  0  for  m=l:M  ,  X  real  ,  (5) 

”  X 

specifies  the  location  of  a  real  SP  in  R^.  This  set  of  M 
simultaneous,  nonlinear  real  equations  must  be  solved  numerically 
for  real  X  e  R^;  if  a  solution  exists,  X  =  X(z)  is  a  function  of 
the  field  point  z  of  interest. 

If  equations  (3)  and  (4)  are  combined,  there  follows 

'  TTTTm  J  exp[A(X,2)]  ,  /g) 

where 

A(X,2)  =  In  T(X,2)  =  In  /j{\)  -  x"^  z  s  x(X)  -  x"^  z  ,  (7) 

and  where  x(X)  is  the  joint  cumulant  generating  function  (CGF) 
corresponding  to  joint  MGF  //(X).  The  Hessian  matrix  (HM)  of 
joint  CGF  X(X)  is  the  MxM  matrix  of  second-order  partial 
derivatives  of  x(X)  with  respect  to  {X^},  m=l:M,  which  are  the  M 
complex  components  of  MD  vector  X.  it  can  be  seen  from  equation 
(7)  that  the  HM  of  A(X,z)  is  identical  to  that  of  X(X),  and  is 
independent  of  field  point  z.  It  is  shown  in  appendix  B  that  the 
HM  of  joint  CGF  X(X)  is  positive  definite  (PD)  for  X  real  and  in 
Therefore,  the  HM  of  A(X,z)  in  equation  (7)  is  also  PD  for  X 
real  and  in  R^.  This  means  that  A{X^,z)  has  positive  curvature 
for  all  Xj.  e  R^;  that  is,  A(X^,2)  has  a  bowl-like  behavior  and 
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More  details  are 


can  have,  at  most,  one  minimum  for  t  R^. 
presented  in  appendix  C. 

From  equation  (7),  since  the  integrand  of  equation  (6)  is 
related  to  A{\^,z)  by  the  monotonically  increasing  transformation 
Y(Xj,,z)  =  exp[A(Xj.,z)  ] ,  the  same  bowl-like  property  carries  over 
to  the  integrand.  That  is,  real  positive  integrand  Y(X^,z)  can 
have,  at  most,  one  minimum  for  X  e  R  .  The  location  of  this 
minimum  (if  it  has  zero  slopes  as  required  by  equation  (5))  is 
the  unique  real  SP  of  Y(X,z)  in  the  ROA(/t/).  This  real  SP  X(z) 
can  be  used  to  connect  the  valleys  of  the  integrand  Y(X,z)  at 
X  *  ±i®  by  means  of  BCs.  Then,  an  SPA  to  p(z)  at  field  point  z 

yv 

can  be  developed  about  this  real  SP  X(z)  in  MGF  space. 

There  may  be  other  complex  SPs  of  integrand  Y(X,z)  for  X  in 
the  ROA(/u);  a  one-dimensional  example  is  furnished  by  PDF  p(z)  = 

A  exp(-z)  [1  +  cos(az)]  for  z  >  0.  For  a  =  1  and  z  =  0.1,  the 

A 

integrand  Y(X,z)  has  a  real  SP  at  X(z)  =  -8.8995  and  a  pair  of 
complex  SPs  at  0.4951  +  i  0.70305,  all  of  which  are  in  the 
ROA(//);  namely,  X^  <  1.  The  AC  of  Y(X,z)  has  another  pair  of 
complex  SPs  at  X(z)  =  1.4547  +  i  0.68369  that  are  outside  the 
ROA(//).  More  details  on  this  example  are  furnished  in  appendix  D. 

For  some  values  of  z,  a  real  SP  of  Y(X,z)  may  not  exist  in 
the  ROA(//).  For  example,  for  one-dimensional  PDF  p(z)  =  exp(-z) 
for  z  >  0,  then  MGF  /y(X)  =  1/(1  -  X)  for  X^,  <  1 .  Then,  if  field 
point  z  is  taken  negative,  z  <  0,  the  real  integrand  Y(X^,z)  is  a 
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strictly  monotonically  increasing  function  of  for  all  X^  in 
(—,1)  and  never  has  zero  slope  there.  This  particular  value  of 
z  IS  unfeasible;  that  is,  z  corresponds  to  a  zero  value  for  the 
PDF  p(z).  The  case  of  z  »  0+  has  an  SP  at  X^.  =  -®.  Figure  1 
illustrates  this  behavior  for  function  A{X^,z). 

But,  even  if  field  point  z  is  a  feasible  point,  the  integrand 
Y(X,z)  may  still  not  have  a  real  SP  in  R^,  even  though  Y{Xj.,z) 
has  a  unique  minimum  for  X^  c  R^.  a  one-dimensional  example  of 
this  situation  is  furnished  by  the  PDF  p(2)  =  a  exp(-z)/(l+z)® 
for  z  >  0.  More  details  on  this  example  are  given  in  appendix  E. 
Thus,  each  candidate  real  solution  of  equation  (5)  must  be 
tracked  as  the  search  for  the  unique  real  SP  in  R^  develops 

in  order  to  ensure  that  the  search  is  not  tending  to  infinity  or 
attempting  to  get  out  of  region  R  . 

fj 


Figure  1.  Logarithm  of  Integrand  for  Exponential  PDF 
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LOCATING  AN  APPROPRIATE  SADDLEPOINT  IN  MOMENT  GENERATING 
SPACE  FOR  EVALUATION  OF  A  NONLINEARITY  AVERAGE 

The  presentation  in  this  section  is  somewhat  similar  to  that 
of  the  previous  section  and  will  be  abbreviated.  Suppose  random 
MD  vector  z  is  subject  to  an  MD  nonlinear  transformation 
g(u)  =  g(Uj^ , . . .  ,Ujj) ,  giving  scalar  output  g(z).  The  average 
of  this  random  quantity  is  given  by 

a  »  E{g(z)}  =  J  du  p(u)  g(u)  ,  (8) 

where  p(u)  is  again  the  joint  MD  PDF  of  RV  z.  Transformation 
g(u)  can  be  a  function  with  additional  parameters,  such  as  the  MD 

.  m 

field  point  z  =  [  z^^  •••  z^]  of  interest,  as  in 
M 

g(u)  =  TT  U(u^  -  zj  ,  (9) 

m=l 

where  U(  )  is  the  unit-step  function,  or  the  power-law  values 

{v„}  in  a  transformation  such  as 
m' 

M  .  V  . 

g(u)  =  TTI%  1  ^m  ^  ®  '  ni=l:M  .  (10) 

Then,  the  average  a  in  equation  (8)  will  be  a  function  of  these 
additional  parameters.  Function  g(u)  is  presumed  to  be  real  and 
nonnegative  for  all  u;  however,  g(u)  need  not  have  finite  volume 
in  u  space,  as  illustrated  by  the  two  examples  above.  Additional 
interpretations  and  examples  are  given  in  reference  2. 
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The  MD  Laplace  transform  of  nonlinearity  g(u)  is  denoted 
y(X)  »  J  du  g(u)  exp(-X^  u)  , 


(11) 


with  a  minus  sign  in  the  exp(  );  this  minus  is  in  contrast  to  the 

Laplace  transform  for  the  joint  MGF  fj(\)  in  equation  (1),  where  a 

plus  sign  was  used.  For  X  purely  real,  that  is,  X.  -  0,  it  is 

presumed  that  integral  (11)  converges  for  X  in  an  MD  region  R 

The  function  y(X^)  is  obviously  positive  real  for  all  X  e  R 

r  Y* 


It  follows  that  if  yCX^,)  exists,  then  y(X)  =  y(X^  +  i  ) 
also  exists  for  all  X^^,  as  seen  from  the  relation 

|y(Xj.  +  i  x.)j  =  |J  du  g(u)  exp[-xj  u  -  i  X^  u]  | 

<  J  du  g(u)  exp[-xj  u]  -  r{\^)  for  all  X.  .  i 


(12) 


This  MD  region  of  complex  X  is  the  region  of  analyticity  of  y(X), 
which  is  denoted  by  ROA(y);  that  is,  ROA(y)  is  the  MD  set  of 
complex  X  values  such  that  Re(X)  e  r^.  Function  y(X)  does  not 
exist  outside  of  ROA(y),  although  its  AC  may.  For  example, 
function  (9)  yields 


T  /  “ 

Y(X)  =  exp(-X  z)/rT(X„)  for  Re(X  ) 

/  m=l  “ 


>  0  ,  m=l;M  . 


That  is,  ROA(y)  here  consists  of  the  right  halves  of  the  M 


complex  planes  {X^j^},  m=l;M. 


(13) 


Equation  (12)  also  indicates  that,  for  a  fixed  real  vector 
e  R^,  the  maximum  value  of  |y(X)|  anywhere  on  the  M  BCs, 

X  =  Xj,  +  i  Xji^,  is  realized  as  all  the  contours  simultaneously 
cross  the  M  real  axes  in  X  space.  Any  perpendicular  departure 
from  the  real  axes,  in  any  or  all  of  the  M  complex  planes 
leads  to  a  decrease  in  |y(X)|  when  X  e  ROA(y). 

Substitution  of  relation  (3)  into  equation  (8)  and  an 
interchange  of  integrals  leads  to  an  alternative  result  for 
average  a;  namely, 

a - f  dX  fj{\)  y(X)  (14) 

(i2ii)”  J 

upon  use  of  equation  (11)  (see  reference  2  for  several  practical 

applications).  However,  this  MD  contour  C  must  now  stay  in  the 

intersection  of  ROAs;  that  is,  X  e  ROA(yu,Y)  =  ROA(yu)  n  ROA(y)  is 

required  in  equation  (14).  Also,  define  the  MD  real  region 

R  =  R  n  R  . 

Mr  /J  r 

The  integrand  in  equation  (14)  is  denoted  as 

Y(X)  *  //(X)  Y(X)  (15) 

and  will  be  a  function  of  any  additional  parameters  that 
nonlinearity  g(u)  depends  on;  for  example,  see  equations  (9)  and 
(10).  If  MD  contour  C  in  equation  (14)  is  taken  as  a  set  of  BCs, 
integrand  Y(X)  will  realize  its  maximum  magnitude  on  C  at 
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the  real  points  where  the  individual  contours  all  cross  the  real 
axes  in  X  space.  Therefore,  the  major  contribution  to  integral 
(14)  is  expected  to  occur  in  the  neighborhood  of  the  points  where 
the  M  contours  cross  the  real  axes  in  X  space.  This  is  the  MD 
point  about  which  to  develop  an  SPA  for  average  a  of  equations 
(8)  and  (14). 

From  equations  (15)  and  (14),  the  logarithm  of  the  integrand 
is  given  by 

A(X)  s  In  Y(X)  =  x(X)  +  In  y(X)  .  (15) 

For  Xj^  =  0,  it  has  been  shown  in  appendix  B  that  the  HM  of  joint 
CGF  X(X)  is  PD  for  e  e  Similarly,  appendix  F 

demonstrates  that  the  HM  of  In  y(X)  is  nonnegative  definite  for 
®  \  ®  Therefore,  the  HM  of  sum  A(X)  in  equation  (16) 

must  be  PD  for  X^.  e  R^^.  This  observation  follows  from  relation 

v^(A+B)v=v^Av+v^Bv,  (17) 

where  v  is  an  arbitrary  real  column  vector,  and  A  and  B  are  real 
square  matrices.  Notice  that  complex  vector  X  must  satisfy  the 
restriction  Re(X)  e  R^^'  '^hile  Im(X)  must  be  zero  for  the  PD 
property  of  the  HM  of  A(X)  to  hold. 

In  the  following  sections,  some  practical  problems  where 
these  issues  have  relevance  will  be  indicated. 
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JOINT  MOMENT  GENERATING  FUNCTION  OF  M  WEIGHTED 
SUNS  OF  N  INDEPENDENT  RANDOM  VARIABLES 

Let  {Xjj}  be  a  set  of  N  independent  real  RVs  with  arbitrary 

first-order  (different)  MGFs  {/;  }  and  corresponding  first-order 

*n 

CGFs  {x„  }.  That  is, 

*n 

(X)  *  exp(X  x^)  ,  X  (X)  =  In  u  (X)  for  n=l:N  ,  (18) 

n  n  *n 

where  the  overbar  denotes  an  ensemble  average.  The  weighted  sums 
of  independent  RVs  of  interest  are  given  by 

N 

y„  '  E  for  m.l!M  ,  (19) 

n=l 

where  M  <  N,  and  MxN  real  matrix  is  arbitrary  except  that 

it  must  have  rank  M. 

The  joint  MGF  of  RVs  {y^^}  is,  upon  use  of  equations  (18)  and 
(19)  and  the  independence  of  RVs 


M 

N 

M 

N  ] 

"y*  ^1  ■ 

}  •  •  m  i 

exp 

n: 

*,iti=l 

\  ^1. 

) 

=  exp 

n: 

^in=l 

n=l  J 

N 

/ 

M 

N 

II 

exp 

*n 

r: 

m=l 

inn  m 

) 

=  n 

n=] 

1  / 

(20) 

re  Mxl 

vector 

X 

tXi 

and 

M 

®  ^  ®mn  ^m  n=l:N  ;  b(X)  *  a'  X  .  (21) 

m=l 
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That  is, 


“  HI  • 


(22) 


The  joint  CGF  of  RVs  using  equation  (18),  is 


N 


yx)  =  i„  ,y(x)  -  c  x,Jb„(x))  . 

There  follows,  by  reference  to  equation  (21), 


(23) 


3X_ 


m 


N 

^  ®mn  ”=1  =  "  • 


(24) 


Finally,  the  M  simultaneous  equations  that  must  be  solved  for  the 
MD  SP  X  are 


N  r  -  ^ 

®mn  '  ^m  ' 


(25) 


where  y  -  [y^  •••  y^^]'  are  the  particular  values  (field  point) 
at  which  the  joint  PDF  of  RVs  {y^j^}  in  equation  (19)  is  to  be 
evaluated.  The  second-order  partial  derivatives  required  to 
obtain  the  SPA  follow  from  equations  (24)  and  (21)  as 


N 


^  Xy(X)  ^  « 

3Xm  3Xjj^  ®mn  ^mn  ^x^yn^^O  m,m=l:M  . 


(26) 


Once  SP  X  is  found  from  equation  (25),  it  can  be  substituted  in 
equation  (26)  and  the  MxM  HM  can  be  evaluated  at  the  SP;  namely. 


'Xy( 


X). 
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EXTENSION  TO  DISTORTED  RANDOM  VARIABLES 


Suppose  that  each  independent  RV  is  the  result  of  a 
different  nonlinear  transformation  according  to 


“  h^(x^)  for  n=l:N  ,  (27) 

where  RV  x^  has  known  first-order  PDF  p„(u).  Then,  the  first- 
— n  ^n 

order  MGF  of  RV  x^  is 

n 

(X)  =  exp(X  Xj^)  =  exp[^X  = 

=  J  du  2jj(u)  exp  ^X  h^{u)j  for  n=l:N  ,  (28) 

from  which  there  follows 


^x  “  J  ^n^^^  exp^X  h^(u)j  for  n=l;N  .  (29) 

If  these  2N  integrals  in  equations  (28)  and  (29)  can  be  evaluated 
in  closed  form,  then  functions 


Xi  (X) 


(X) 

— .  for  n=l:N 
(X) 

*n 


(30) 


are  available  for  use  in  the  determination  of  the  SP  required  in 
equation  (25).  Otherwise,  numerical  integration  could  be  used  on 
these  one-dimensional  integrals. 
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PROBABILITY  DENSITY  FUNCTION  OF  THE  MAXIMUM  RANDOM  VARIABLE 


T 

For  M  =  2,  let  scalar  RV  y  =  niax(  ,  Z2 ) .  Let  u  =  U2  ]  , 

T 

Y  =  [y  y]  f  and  real  nonlinearities 

g^(u)  =  6(y  -  u^)  U(y  -  U2)  ,  92^“^  “  “  ^1^  * 

Then,  total  nonlinearity  g(u)  =  g2(n)  +  92^^)  is  real  and 
nonnegative.  The  average  corresponding  to  total  nonlinearity 
g(u)  is 

y  y 

a  =  J  du  g(u)  P2(n)  =  J  du2  P2(y,U2)  +  J  du^  p^(Uj^,y)  ,  (32) 

^00  —00 

which  is  the  PDF  Py(y)  of  scalar  RV  y  =  max(Zj^,Z2)  at  argument  y. 


The  corresponding  gamma  functions  to  gj^(u)  and  92^^)  are 
T 

Yj^(X)  =  <  0  ,  arbitrary;  (33) 

2 


Y2(X)  =  Re(Xj^)  <  0  ,  X2  arbitrary. 


(34) 


The  total  gamma  function  corresponding  to  nonlinearity  g(u)  is 

y(X)  =  Yi(X)  +  Y2(X)  for  Re(Xj^)  <  0  ^  Re(X2)'  <  0  .  (35) 

However,  notice  that  combining  the  two  individual  gamma  functions 
together  has  forced  a  severe  restriction  on  the  allowed  region  in 
X-space  to  investigate  for  an  SP;  namely. 
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^  *  rt]  .  (36) 

where  component  contours  and  C2  of  must  both  pass  to  the 
left  of  their  origins.  For  y  >  0,  the  restrictions  to  the 
left-half  X  planes  keep  the  real  part  of  the  argument  of  exp(  ) 
positive;  therefore,  exp(  )  is  unable  to  track  the  rapidly 
decreasing  upper  tail  of  the  PDF  of  y.  The  SP  components 
approach  0-  in  all  M  =  2  dimensions  as  threshold  y  increases,  but 
this  is  not  sufficient  to  allow  development  of  a  rapid  decrease 
of  the  SPA  on  the  upper  tail  of  the  PDF.  See  figure  10  of 
reference  2  for  an  example  of  this  limitation. 


On  the  other  hand,  for  nonlinearity  gj^(u)  alone,  the 
contribution  to  average  a  is 


1 

(i2n)^ 


C 


dX 

1) 


exp(-X^y-X2y)/(-X2) 


(37) 


which  requires  only  that  the  X2  contour  pass  to  the  left  of  its 
origin.  Then,  the  Xj^  component  of  the  SP  can  move  into  its 
right-half  plane  and  yield  exponential  decay  of  a^^  with  y,  even 
when  threshold  y  >  0.  Conversely,  for  average  a-,  the  X- 
component  of  the  SP  is  free  to  move  as  needed.  in  this  manner, 
much  better  SPAs  are  obtained,  as  well  as  smaller  correction 
terms.  The  tradeoff  is  that  M  =  2  SPs  must  be  located,  instead 
of  a  single  SP  for  equation  (36). 
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In  M  dimensions,  the  m-th  nonlinearity  required  to  yield  the 
PDF  of  scalar  RV  y  =  max(Zj^, _ ,z^)  is 

M 

gjj^(u)  =  8(y  -  J~J  u(y  -  u^)  for  m=l;M  ,  (38) 

n=l 

n?^m 

for  which  the  corresponding  gamma  function  is 


m 


f  M  N  /  M 

=  exp  -y  IZl^n  /  for  m=l;M 

''  n=l  n=l 


n^^m 


(39) 


provided  that  Re(X  )  <  0  for  n=l:M,  n/m,  while  X  is  arbitrary. 

II  ni 

Therefore,  the  m~th  average 


m 


( i2n) 


M 


dX 

m) 


//^(X) 


exp 


n^m 


(40) 


allows  for  the  m-th  component  contour  to  move  into  the 

right-half  X^^^-plane  (for  each  m  value  in  [1,M])  and  thereby 
realize  significant  exponential  decay  for  a^j^,  even  if  threshold 
y  >  0.  On  the  other  hand,  an  attempt  to  minimize  the 
computational  effort  and  use  total  gamma  function 


TO')  -  C  T„(X)  , 

in=i 


(41) 


as  in  equation  (36)  for  M  =  2,  would  require  that  Re(Xjj^)  <  0  for 
ni*l;M,  and  thereby  severely  limit  the  capability  of  the 
resulting  SPA,  especially  for  y  >  0. 


19  (20  blank) 


JOINT  STATISTICS  FOR  NORMALIZATION 


Random  vector  x  =  [x^^  •••  x^j] '  has  arbitrary  joint  PDF  and 
joint  MGF  The  individual  RVs  {Xj^^}  can  be  positive  or 

negative,  and  can  be  statistically  dependent  on  each  other. 


GENERAL  NORMAL I ZER 


Scalar  RV  y  may  be  statistically  dependent  on  random  vector 
x;  however,  y  is  positive  (for  the  time  being).  The  first-order 
PDF  of  y  is  Py.  The  normalizer  of  interest  forms  the  M  ratios 


z  =  —  for  m=l:M  . 
m  y 


(42) 


The  joint  cumulative  distribution  function  (CDF)  of  random  vector 
z  =  [z^  •••  JBjjl'  using  y  >  0  is 

c^(Zi,  —  ,Zj^)  =  Prob(Zj^  <  z^,  —  <  z^) 

=  Prob(Xj^  <  z^  y,...,Xjj  <  z^^  y)  = 

»  ^ly  ^M^ 

=  J  dy  J  dxi  •••  J  dXj^  p^y(x^,...,Xj^,y)  - 

Q  .CD 

^ly  ^My 

=  J  dy  Py(y)  J  dxj^  •••  J  dx^^  p^(x^, . . .  ,Xjj|y  =  y)  ,  (43) 

Q  .00  .00 

where  p  is  the  joint  PDF  of  x  and  y,  and  p  (x|y  =  y)  is  the 
xy  X 

conditional  joint  PDF  of  RV  x  at  argument  x,  given  that  scalar  RV 
y  =  y. 
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The  joint  PDF  of  normalizer  vector  output  z  follows  from 
equation  (43),  upon  M  partial  differentiations,  as 

CO 

=  J  dy  y  py(y)  p^(z^  . . 2^  y|y  =  y,  ^  (44^ 

0 

00 

'  M 

~  y  Pxyt^i  y>  —  y»y)  .  (45) 

0 

If  joint  PDF  p^y  can  be  determined,  a  one-dimensional  integral 
must  be  evaluated  in  order  to  obtain  joint  PDF  p 

z 

If  ^  the  RVS  {Xjjj}  are  statistically  independent  of  scaling 

y  used  in  normalizer  ratios  (42),  then  equation  (44)  simplifies 
to 

00 

Pz<^1'-*-'2m>  =  j  dy  y”  Py(y)  p^(z^  y)  .  (46) 

0 

This  result  still  requires  a  one-dimensional  integral  for  its 
evaluation . 

If  RV  y  can  take  on  both  positive  and  negative  values,  the 
three  integrals  in  equations  (44)  through  (46)  are  generalized  to 
the  extent  that  the  integrals  over  y  are  conducted  from  -«  to  +», 
and  the  term  y«  is  replaced  by  |y|«.  Determination  of  the 
conditional  PDF  in  the  integrand  of  equation  (44)  can  be  very 
involved;  it  may  be  safer  and  easier  to  use  the  form  involving 
the  joint  PDF  in  equation  (45).  in  general,  the  analytic 
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difficulty  of  determining  joint  PDF  p  prompts  us  to  consider 

xy 

specific  forms  of  normalizers,  and  attempt  to  find  the  joint  MGFs 
of  their  outputs. 


A  SPECIFIC  NORMALIZER 

Let  denominator  scaling  y  in  normalizer  ratios  (42)  be  formed 
by  linear  weighting  of  random  vector  x  according  to 

M 

y  C  (47) 

where  scalar  RV  w  is  independent  of  x.  For  example,  w  might 
involve  averaging  additional  RVs  {Xj^^},  m  >  M,  which  are 
independent  of  {x^},  m=l!M.  The  MGF  of  scalar  RV  w  is  fj^. 

Observe  that  numerators  {Xjjj}  are  statistically  dependent  on 
denominator  y  in  normalizers  (42)  and  (47);  the  one  exception  is 
when  a  =  0  in  equation  (47). 

It  is  presumed  that  scalar  RV  y  is  always  positive.  Then, 
upon  substitution  of  equation  (47)  into  equation  (42),  the  joint 
CDF  of  normalizer  output  z  is 

. "«’  ■  . . - 

■  y  -  >  0 . y  -  *„  >  0)  - 

-  Prob(Vj  >  0 . >  0)  ,  (48, 

where  RVs 
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The  joint  MGF  of  random  vector  v  is 
A/y(X)  =  exp(  \'v)  *  exp(X'Ax  +  X'zw)  = 

=  exp(x'  b(X))  exp(  X'zw)  =  f/^{h{X))  /j^{z'\)  ,  (52) 

using  the  independence  of  random  vector  x  and  random  scalar  w, 
and  defining  Mxl  vector 

b(X)  =  A'X  ;  b^{X)  =  YU  \  for  m=l;M  .  (53) 

ni=l  ~  — 

It  is  also  presumed  above  that  b(X)  e  ROA(//^)  and  that 

z'X  e  ROA(//^).  Finally,  the  joint  CDF  of  normalizer  output  z 

follows  from  equations  (48)  and  (52)  as  (see  reference  2) 
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(54) 


In  the  special  case  where  RVs  are  independent  of  each 

other,  but  not  necessarily  identically  distributed,  the  joint  MGF 
required  in  equation  (54)  becomes 


(55) 


in  terms  of  the  individual  MGFs  {//^  }  of  RVs  Notice  that 

in 

the  statistical  dependence  of  all  the  RVs  n=l;N,  on 

denominator  RV  y  in  equations  (42)  and  (47)  is  fully  accounted 
for  in  equations  (54)  and  (55),  basically  through  MxM  matrix  A  in 
equation  (51)  and  Mxl  vector  b(X)  in  equation  (53). 


For  general  RVs  if  the  additional  RV  w  in  normalizer 

denominator  y  in  equation  (47)  is  zero,  then  yu^(u)  =  1  for  all  u; 
this  simplifies  the  evaluation  of  equation  (54).  However, 
observe  that  A  =  A(z,a),  leading  to 

b(X)  =  A(z,a)'  X  =  b(z,a,X)  ,  (56) 

through  equation  (53).  On  the  other  hand,  if  w  is  a  nonzero 
constant,  then  its  MGF  is  simply  yu^(u)  =  exp(wu),  thereby  leading 
to  the  term  exp(wz'X)  in  equation  (54)  for 
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Alternatively,  if  scale  factor  a  =  0  in  equation  (47),  then 
matrix  A  in  equation  (51)  is  -i^,  and  equation  (53)  yields 
b(X)  =  -X.  Then,  equation  (54)  yields  joint  CDF 


2'^)/rT(x„) 


(57) 


If  RVs  (x^)  are  independent  of  each  other,  there  follows 


^^(-X)  .  ^  (-X^)  . 

m=l  m 


(58) 


ALTERNATIVE  PDF  APPROACH 


In  equation  (49),  define  RVs 


(59) 


Since  this  is  a  one-to-one  transformation,  the  joint  PDF 
be  immediately  found  from  joint  PDF  p^.  Then,  for  w  >  0, 
equations  (48)  and  (49)  take  the  form 

c^(z)  =  Prob(u^  >  . . u„  >  -z„  w)  = 

00 

=  J  dw  p^(w)  e^(-z^  w,  ...,-Zjj  W)  , 


PDF  p^  can 


(60) 


using  the  independence  of  vector  u  and  RV  w,  as  well  as  the  joint 
exceedance  distribution  function  e„.  The  corresponding  Joint  PDF 
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00 


p^(z)  = 


dw  w' 


M 


P„(w) 


P  (“Z, 
1 


w. 


,-z 


M 


w) 


0 

which  agrees  with  equation  (46)  under  a  re-identification  of 
variables . 


(61) 
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IMPROVING  THE  ACCURACY  OF  THE  H-DINENSIONAL 
SADDLEPOINT  APPROXIMATION  BY  MEANS  OF  GAUSS  QUADRATURE 


Let  X  = •••  Xjj]'  be  an  MD  RV  with  joint  MGF 
/Ujj(X)  =  exp(  X'x) ,  where  vector  X  =  [  X^  •••  Xj^]'.  Then,  the  joint 
PDF  of  X  at  MD  field  point  x  =  [x^  •••  is  given  by  the  M-th 

order  integral 

^  f  exp(-X'x)  /u  (X)  ,  (62) 

where  BC  C  is  parallel  to  the  imaginary  axis  in  each  of  the  M 
dimensions.  The  joint  CGF  of  x  is  Xjj(X)  =  In  •  The 

MD  SP  of  the  integrand  of  equation  (62)  is  that  real  point 
X  *  X(x)  nearest  the  origin  where  all  M  partial  derivatives 
satisfy 


3X  (X) 

for  m.l!M  .  (63) 

"  X 


When  the  contour  C  is  moved  in  M  dimensions  so  as  to  go  through 
the  SP,  and  the  change  of  variable 

A. 

X  =  X  +  i  t  ,  t  =  [t^  ••  •  tjjl'  (64) 

is  made,  equation  (62)  becomes 

Px(x)  =  (2ji)“”  exp(-X'x)  J  dt  exp(-it'x)  //^(X+it)  ,  (65) 

where  the  new  contour  passes  through  the  SP  of  the  integrand  of 
equation  (65),  which  is  now  at  t  =  0. 
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The  logarithm  of  the  integrand  of  equation  (65)  can  be 
expanded  in  a  power  series  about  the  origin  according  to 


ln{exp(-it'x)  //^(X+it)}  -  -it'x  + 


Xjj(X+it)  = 


=  -it'x  +  x„(X)  +  it'  VX^{X)  -  4 


2  t'  H  (X)  t  , 


(66) 


where  MxM  matrix 


„  ,,,  _ 

*  3X - 3jr  '  , 

m  m. 


(67) 


is  symmetric  im  m  and  m  for  all  X.  Thus,  using  equations  (66) 
and  (63),  the  integrand  of  equation  (65)  can  be  approximated  as 


exp(  it'x)  //^(X+it)  :  exp(x^(X)  -  |  f  Hjj(X)  t] 


(68) 


for  small  jt|.  if  this  approximation  is  now  extrapolated  to  all 
t  and  substituted  in  equation  (65),  there  follows  the  usual  SP 
(or  tilted  Edgeworth)  approximation  in  M  dimensions: 


Px(x) 


~  ®^p(Xx(X)-X'x] 


( 2n) 


J  dt  exp(-  I  t'  H^(X)  t]  = 


exp[Xx(X)-X'x] 

det(H^(X)]^  '*  "  =  • 


(69) 


The  extrapolation  of  approximation  (66)  to  all  t  was  done  foi 
purposes  of  analytically  evaluating  the  MD  integral  in  equation 
(65).  If,  instead,  the  MD  change  of  variable 
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(s  is  Mxl) 


(70) 


t  =  H^(X)  s 


is  made  in  equation  (65),  there  follows  the  exact  PDF  form 


p^(x) 


exp  1 

f  ^ 

X„(X) 

^  * 

,-i-x) 

(2ii)”/^  det 

H,(X) 
^  * 

(71) 


where  the  multiplicative  M-fold  integral  is 


(2 


]  as  exp(-i  s>  xl  . 


(72) 


The  integrand  of  this  latter  integral  behaves  as  exp(-s's/2)  near 
the  SP  at  the  origin  of  s  space;  this  leads  to  the  earlier 
approximation  above,  namely,  1=1. 


2 

For  one-dimensional  integrals  in  x  that  behave  as  exp(-a  x  ) 
near  the  origin,  the  use  of  Gauss  quadrature  (reference  3,  page 
924)  suggests  itself  as  a  good  candidate  for  high  accuracy  with 
few  integrand  evaluations.  This  quadrature  approach  can  be 
readily  extended  to  M  dimensions  by  simply  repeating  the 
one-dimensional  rule  for  each  of  the  M  dimensions  under 
investigation.  This  is  consistent  with  the  fact  that  integral 
(72)  has  been  transformed,  and  for  which  the  integrand  behaves 
similarly  in  all  dimensions,  namely  as  exp(-Sjj^/2)  in  the  m-th 
dimension,  at  least  near  the  SP  at  the  origin  s  =  0.  The 
numerical  question  to  answer  in  practice  is:  how  close  is 
integral  I  in  equation  (72)  to  the  value  1? 
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The  basic  problem  with  the  quadrature  approach  in  high 
numbers  of  dimensions  is  that  if  S  samples  are  taken  per 
dimension  of  the  integrand,  the  total  number  of  samples  required 
in  M  dimensions  is  s”.  For  example,  just  S  =  2  samples  per 
dimension  requires,  in  M  =  16  dimensions,  66,000  evaluations  of  a 
complex  integrand  which  is  a  function  of  M  variables.  This  can 
be  a  very  time-consuming  task  and  require  considerable  storage, 
depending  on  how  large  M  gets. 


Define  the  residual  integrand  of  equation  (72)  as  the  actual 
integrand  multiplied  by  exp(s's/2);  this  residual  integrand  is 
presumed  to  be  1  for  all  s  under  the  extrapolation  assumption 
leading  to  equation  (69).  On  the  other  hand.  Gauss  quadrature, 
using  S  samples  in  a  particular  dimension,  is  effectively  fitting 
the  residual  integrand  with  a  polynomial  of  order  2S-1  in  that 
variable,  which  affords  a  rather  high-order  fit  for  few  samples, 
in  M  dimensions,  it  means  that  the  Gauss  quadrature  procedure  is 
for  residual  integrands  of  the  form 


for  any  constants  {a(P2  f  • .  »Pjj) } ,  where  integer  p^  <  2S-1  for  all 
m-l:M.  This  potentially  high  order  of  accuracy  makes  the  use  of 
Gauss  quadrature  on  equation  (72)  a  strong  candidate  for 
increased  accuracy  of  the  SPA  if  the  large  number  of  integrand 
evaluations  can  be  tolerated. 
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FIRST-ORDER  CORRECTION  TERM  TO  SADDLEPOINT  APPROXIMATION 

The  first-order  (FO)  correction  term  to  the  SPA  (or  tilted 
Edgeworth  approximation),  that  is,  the  0(1/N)  term,  is  given  in 
reference  4,  page  180.  To  develop  it  fully  for  these  purposes, 
the  following  definitions  are  employed  here.  For  RV  x  with  MD 
CGF  Xjj/  and  MD  field  point  x  of  interest  with  corresponding  SP 
X  =  X(x),  let  j,k,X,m  =  1;M,  and  define 

3^Xx(X) 

^Xm  “  axTTxI  .  ' 

X  m 
3^Xx(X) 

"  ax.  3X,  .  ' 

K  A  ^  X 

3^Xx(X) 

^jkxm  ”  ax.  ax.  ax,  ax^  >.  • 

3  K  A  nt 

Also,  define  the  two  MxM  matrices 

H  =  tX^l  ,  T  -  IT^l  -  H-1  .  (77) 

H  is  the  HM  of  the  MD  CGF  Xj^  of  RV  x,  evaluated  at  the  SP  X. 

Then,  the  FO-corrected  joint  PDF  at  MD  field  point  x  is 
p^(x)  =  H  +  +  Cj^  +  Cjjj]  »  (78) 

where 


(74) 


(75) 


(76) 
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expfx,(X)  -  X'  xl 

Pn(x)  »  -r-^ - 1- 

((2n)“  det(H)]^ 

IS  the  usual  { zeroth-order )  SPA,  and  parameters 

®  jSi  ’’jk  ’’*1"  ' 


-  -  J  C  C  X 


kXm  kXm 


kXm  ^kXm  "^kX  "^mk  ’’xm  ' 


'3b  "  12  ^klm  Xv,_  T..  T,-  T 

kXm  kXm  XX  i 


(82) 


terms  (c^,  ‘^Sa'  0(1/n)  if  n  independent  RVs 

have  been  added  te  form  RV  x.  However,  the  results  above  hold 
for  any  RV  x,  no  matter  how  formed.  The  two  component  terms, 
Cja  and  cannot  be  combined  together.  In  general. 


APPLICATION  TO  WEIGHTED  SUMS  OF  RANDOM  VARIABLES 


Let  N-dlmensional  RV  w  .  Iw^  ...  be  a  collection  of  N 

independent  RVs,  with  possibly  different  statistics  for  each  RV 
n=l:N.  in  particular,  let  individual  RV  w„  have  PO  MGP 
and  FO  CGF  for  n.l,N.  The  MD  RV  x  of  Interest  is  formed  "" 
according  to  linear  transformation 


X  =  a  w  ,  a  =  fa  1  . 

L  mnj 

Arbitrary  array  a  is  MxN,  M  <  n,  and  of  rank  M.  Then,  Rv  x  i 


(83) 
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composed  of  M  dependent  RVs  with  joint  MGF 


|/3j(X)  =  exp(X'  x)  =  exp(X'  a  w)  =  exp(b(X)'  w)  ,  (84) 

where  Nxl  vector 

b(X)  =  a'  X  =  [b^(X)  •••  bj^(X)]'  ,  (85) 

with  elements 

M 

=  rZ  a^n  ^m  '  (86) 

m=i 

Notice  that 
3b^(X) 

"ax -  =  ®mn  •  (87) 

m 

Using  the  independence  of  RVs  {w^},  there  follows  for  the 
joint  MGF  of  X  from  equations  (84)  through  (86), 

■  n  ''n]  -  TJ  ■  (88) 

The  joint  CGF  of  RV  x  is  then 

X,(^)  -  C  •  (89) 

There  follows,  upon  use  of  equation  (87), 

9Xx(X)  N 

S  '  (98) 
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3  X3^(X)  N 


X=1:M  ,  m=l:M  .  (91) 

Therefore,  the  M  simultaneous  equations  that  must  be  solved  for 
the  SP  X  -  X(x)  at  MD  field  point  x  are 


^mn  “  ^m  in=l:M 


(92) 


In  order  to  develop  explicitly  the  component  terms  in 
equations  (80)  through  (82),  define,  for  n=l:N, 

’'2n  -  .  X3„  -  X;^(b„(X))  ,  X4„  -  X;”(b„(M)  .  (93) 

These  3N  quantities  only  need  to  be  computed  at  the  SP  X.  The 
use  of  equations  (93)  and  (74)  in  equation  (91)  yields 


^Xm  ^Xn  ®mn  ^2n  X=1:M  ,  m=l;M  . 


(94) 


in  a  similar  fashion,  there  follows  from  equations  (75),  (76), 
and  (93),  for  j,k,X,m  -  IsM, 


Xk*m  “  IZI 


®kn  ®Xn  ®mn  ^3 


(95) 


IN 

^jkXm  ®jn  ®kn  ®Xn  ®mn  ^4n  ‘ 


(96) 


The  three  component  terms  in  equations  (80)  -  (82)  can  now  be 
determined.  From  equations  (80)  and  (96),  there  follows 
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8  ^ — '  ®jn  ®kn  ®Xn  ®mn  ^4n  ^jk 

3  k  Xm  n 


8  ^4n 

n=i 


[jS  "Jn  ^jk  *knj 


Define  the  NxN  matrix  J  =  with  elements 


M 

=  HH  a„  for  n,n  =  1:N 


Then, 

j  =  a'  T  a  , 

and  equation  (97)  can  be  expressed  as 

1  N  2 

^4*8  ^ i  ^4n  "^nn  ’ 
n=l 

This  result  uses  only  the  diagonal  of  NxN  matrix  J. 


Also,  from  equations  (81)  and  (95),  there  follows 

^3a  ”  ”  F  ^-s-*  ^ — '  ^kn  ®Xn  ®mn  ^3n  ^kn  ^Xn  ^mn  ^3n 

kXm  kXm  n  n  —  —  — 

1  ” 

^kX  "^mk  "^Xm  ~  ~  8  ^3n  ^3n  “^nn  "^nn  "^nn 

—  —  nn=i  —  —  — 

Define  Nxl  vector  j  =  [jj^  •••  jjjl^  with  elements 


j 

Then, 


n 


-  ^nn  "-1=”  ' 

equation  (101)  can  be  expressed  compactly  as 


(97) 


(98) 


(99) 


(100) 


(101) 


(102) 
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1 

8 


(103) 


c 


3a 


N 


c 

nn=l 


'^nn 


'n 


J  j  . 


The  final  component  term  is,  from  equations  (82)  and  (95), 


'3b 


^3n  C  ® 


kn  ®Xn  ^mn  ^3n  * 


"^kk  ^XX  ^mm 


N 


1 

-  T7  e:  X 

nn=l 


12  ^3n  ^3n  "^nn  • 


(104) 


The  total  correction  term  of  0(1/N)  required  for  p^(x)  in 
equation  (78)  is  given  by  the  sum  of  components  (100),  (103),  and 
(104).  The  inputs  required  here  are  MxN  matrix  a,  individual  FO 
CGFs  {x^^}  of  RVs  {w^}  for  n=l;N,  and  field  point 
X  =  [x^  ...  x„]'. 


The  MATLAB  code  for  accomplishing  all  of  these  matrix 
manipulations  is  very  compact  once  the  Nxl  chi-vectors  in 
equation  (93)  are  available: 


b*a' *lambda ; 

H»a.*repmat(chi2' ,M,l)*a" ; 

J-a  VH*a,. 

Jd-diag( J ) ; 
j=chi3 . *Jd; 
c4=chi4'*(  Jd.'‘2)/8,. 
c3a=-j  '*j*j/8,. 
c3b=-chi3'*(  J.^3)*chi3/12,. 
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EXACT  H-DIHENSIONAL  PROBABILITY  DENSITY 
FUNCTIONS  FOR  SPECIAL  TRANSFORMATIONS 

Let  X  be  a  column  vector  of  M  real  RVs  with  joint  PDF  p^. 
Also,  let  the  real,  nonrandom  MxM  matrix  A  be  of  rank  M.  Then, 
the  joint  PDF  of  the  column  vector  of  M  real  RVs  y  =  A  x  can 
be  determined  immediately  as 

Py(y)  =  y) ^ |det(A)  I  ,  (105) 

where  y  =  [y^^  ...  y^j]'  is  the  nonrandom  field  point  of  interest. 
It  is  desired  to  extend  this  exact  result  for  the  PDF  of  RV  y  to 
matrices  A  of  size  MxN  of  rank  M,  where  N  may  be  (much)  larger 
than  M.  This  extension  is  not  achievable  for  all  possible  A 
matrices,  but  can  be  accomplished  to  yield  exact  PDFs  of  RV  y  for 
linear  transformations  A  of  a  particular  form. 

Let  V( 1 ) , . . . , V(M)  be  arbitrary  Mxl  real  column  vectors  that 
are  linearly  independent  of  each  other.  Form  the  MxN  matrix 

A  s  [d^  V(m(l))  V(m(2))  •••  d^j  V(m(N))]  ,  (106) 

where  N  >  M  and  d^,  d2».../  d^j  are  N  arbitrary  real  scalars, 
while  m(l),  m(2),...,  m(N)  are  N  arbitrary  integers  in  the  range 
[1,M].  To  make  the  rank  of  A  equal  to  M,  every  vector  V(l) 
through  V(M)  must  appear  at  least  once  in  matrix  A  in  equation 
(106).  Alternatively  stated,  the  set  of  N  integers  m(l),  m(2), 
...,  m(N)  must  include  every  integer  in  the  range  [1,M]  at  least 
once . 
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Let  X  =  •••  X  ]',  where  the  N  scalar  RVs  {x  }  are 

statistically  independent  of  each  other  but  are  not  necessarily 
identically  distributed.  Consider  the  Mxl  RV  obtained  by  the 
linear  transformation 


N 

^  '  S  v(m(n))  Xn  -  ^^1  •••  ^  (107) 

where  explicit  use  of  equation  (106)  has  been  made.  Now, 
let  denote  the  set  of  n  values  where  m(n)  =  1; 

let  Sjj  denote  the  set  of  n  values  where  m(n)  =  M.  (108) 

Then,  equation  (107)  can  be  expressed  as 


^  '^n  •••  +  V(M)  j;  d  X 

where  the  M  RVs  {w^^}  are  defined  as 
nes_ 


m=l 


V(m) 


w_ 


m 


(109) 


(110) 


The  M  scalar  RVs  are  statistically  independent  of  each 

other  because  sets  S^,...,Sj^,  defined  in  equation  (108),  are 
disjoint  and  non-empty;  therefore,  no  common  RVs  of  set  {x  >  can 
appear  in  two  different  RVs  of  set  .  The  FO  characteristic 

function  (CP)  of  the  m-th  RV  w^  follows  immediately  from  equation 
(110)  as 


f  (5) 

m 


TT  (d„c) 

neS  n 
m 


for  m=l:M  , 


(111) 
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upon  use  of  the  independence  of  the  M  RVs  {w  }  defined  in 

ni*' 

equation  (110). 
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In  order  to  utilize  equation  (116),  it  is  necessary  that  each 
FO  PDF  capable  of  evaluation  at  arbitrary  argument  w; 

but  this  quantity  is  available  by  means  of  a  one-dimensional 
Fourier  transform  of  equation  (111): 

Pw  '  2rt  J  exp(-iw^)  f  (5)  for  m=l:M  .  (117) 

®  m 

For  a  given  or  specified  set  of  values  y  =  (y^,...,yj^)  in 
equation  (116),  argument  =  b^^  y^  +  . . .  +  y„  for  PDF  p^ 
can  be  computed,  and  single  integral  (117)  can  be  evaluated  for 
this  particular  value.  This  basic  integral  must  be  redone  for 
each  m  in  the  range  m=l:M  for  use  in  equation  (116).  An  exact 

value  for  the  PDF  Py(y)  of  RV  y  in  equation  (107)  can  be  obtained 
in  this  fashion. 


SPECIAL  CASE 


Let  the  N  integers  m(l),  ...,  m(N)  in  equation  (106)  consist 
of  K  Is  followed  by  K  2s  ...  followed  by  K  Ms,  for  a  total  of 

N  -  KM  integers.  Then,  set  in  equation  (108)  is  the  set  of 
integers 


=  {l+(m-l)K, . 

Use  of  this  relation 

mK 

m  n=l+mK-K 


. ,  mK}  for  m=l :M  . 
in  equation  (111)  yields  CF 

m=l:M  . 

*n  " 


(118) 


(119) 
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Furthermore,  specialize  to  the  case  of  independent 

identically  distributed  RVs  {x^^}  with  an  exponential  PDF  with 

unit  mean;  that  is,  their  common  CF  is  (^)  *  (1  -  i^)  ^  for 

n 

n=l:N.  Then,  equation  (119)  yields  the  explicit  CF  relation 


it) 


mK 

TT 

U=l+mK-K 


for  m=l:M  . 


(120) 


This  result  must  now  be  subjected  to  integral  (117)  in  order  to 
determine  the  FO  PDFs  of  RVs  . 


As  a  further  specialization,  consider  sequence  d^^  =  1  for 

n=l;N.  Then,  equation  (120)  yields  CF  f  ( ^)  =  ( 1  -  i ^)  ^  for 

m 

all  m=l:M,  and  equation  (117)  yields  closed-form  PDFs 


P„  (w)  - 
m 


K— 1 

w  .  exp(-w) 
(K-D! 


U(w)  for  m»l:M 


(121) 


where  U  is  the  unit-step  function.  This  result  can  be  applied  in 
equation  (116)  for  the  joint  PDF  Py(y)  upon  use  of  B  =  v“^. 


Alternatively,  consider  the  complete  scalar  sequence  {d^^}  to 

be  composed  of  M  groups  of  K  scalars  each.  Let  each  group 

contain  J  (=  K/2)  +ls  and  J  -Is,  K  even.  Then,  equation 

(120)  yields  CF  f  (^)  =  (1  +  52)-K/2  =  +  ^2^-0  ^^Ll  m-l;M, 

m 

and  equation  (117)  yields  FO  PDFs 


Pw 

m 


^  r2J-2-j'i 


for  m=l:M 


(122) 


for  all  w.  Again,  direct  application  in  equation  (116)  gives 
exact  results  for  joint  PDF  Py(y)  at  any  argument  y. 
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SUMMARY 


The  dominant  saddlepoint  in  the  M-dimensional  region  of 
analyticity  of  the  original  moment  generating  function  is  located 
on  the  real  axes  of  X  space,  if  it  exists  at  all.  The  numerical 
search  for  this  saddlepoint  is  alleviated  by  the  fact  that  the 
real  integrand  of  interest  has  positive  curvature  (a  bowl-like 
behavior)  in  the  M  real  dimensions  with  a  single  minimum  in  the 
region  of  analyticity.  The  first-order  correction  term  to  the 
standard  saddlepoint  approximation  requires  the  calculation  of 
numerous  fourth-order  partial  derivatives  of  the  joint  cumulant 
generating  function  of  the  random  variables  of  interest;  however, 
these  calculations  can  be  limited  to  the  saddlepoint  location 
alone . 
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APPENDIX  A  —  BRANCH  POINT  EXAMPLE 


Consider  the  one-dimensional  PDF  and  MGF  pair: 

“  >  0  .  v  >  -1  -• 


fj{\)  =  (1  -  X)  ^  ^  for  <  1 

/u(Xj.)  is  positive  real  for  X^  <  1. 
at  X  =  1  when  v  is  not  an  integer, 
the  PDF  expression 


(A-1) 

The  AC  of  //(X)  has  a  BP  at 
Use  this  AC  of  fj{\)  to  get 


“  Tlii  J  exp(-uX)/(l-X)^‘^^  ,  (A-2) 

C 

where  BC  C  passes  to  the  left  of  X  =  1. 


Let  u  >  0.  Also,  let  —1  <  v  <  0  for  now.  Wrap  contour  C 
around  the  BP  at  X  =  1,  making  a  keyhole  contour  centered  on  the 
posi tive— real  X— axis;  this  equivalent  contour  is  the  steepest 
descent  contour.  The  infinitely  small  circular  contour  around 
X  *  1  yields  a  zero  contribution  in  the  limit  because  v  <  0.  Let 
^  1  exp(i<^);  1  -  X  =  -r  exp(i(|>)  =  r  exp( i<|)-in) . 

When  ♦  =  11,  1-X=r>0  and  (1-X)^^^  is  positive  real,  as 
required.  On  the  upper  line  integral,  ♦  =  0,  giving  (l-X)'^'*’^  = 
r  exp t  — i ji ( v+l )  ] ;  on  the  lower  line  integral,  ♦  =  2ii,  giving 

®  r^'*’^  exp[ i n( v+1 )  ] .  Substitution  into  equation  (A-2) 
yields,  for  u  >  0, 


A-1 


on  the  other  hand,  the  SPA  uses  the  SP  of  the  integrand  of 
equation  (A-2);  namely,  Y(X,u)  =  exp(-uX)/(  at  location 
X(u)  »  1  -  (v+l)/u  for  u  >  0.  The  resulting  SPA  is 


rN  _  exD(v+l)  V 

°  “  (2n)^  for  u  >  0  ,  V  >  -1  .  (a-4) 


This  form  is  identical  to  PDF  p(u)  in  equations  (A-1)  and  (A-3) 
except  for  a  scale  factor.  The  ratio  Pq(u)/p(u)  is  independent 
of  u  and  approaches  1+  as  v  increases;  the  ratio  is 
e//(2n)  -  1.0844  at  v  =  0,  and  is  eV(4/ii)  *  1.0422  at  v  »  1 . 
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APPENDIX  B  -  HESSIAN  MATRIX  OF  JOINT  CUHULANT  GENERATING 

FUNCTION  X(X)  IS  POSITIVE  DEFINITE  FOR  X  REAL  AND  IN  R 

f! 

Vector  field  point  z  =  [  z^^  •••  is  real  and  MD.  Let  p(z) 

be  the  joint  PDF  of  MD  RV  z  at  field  point  z;  thus^  p(z)  is  real, 
nonnegative,  and  has  unit  volume  in  MD  z  space. 


T 

Vector  X  =  •••  X^^]  is  complex  and  MD.  The  joint  MGF 

//(X)  corresponding  to  joint  PDF  p(z)  is  given  by  Laplace 
transform  and  expectation 

//(X)  =  J  dz  exp(X^  z)  p{z)  =  Efexptx"^  z)}  .  (B-1) 

Henceforth,  it  is  presumed  that  X  e  ROA(/y);  that  is,  X^,  e  R^. 
Define 

g 

^(X)  =  yWjjj(X)  for  m=l:M  .  (B-2) 

in 

Then,  the  joint  CGF  X(X)  of  RV  z  satisfies  the  relations 


X(X)  =  In  , 


ax 


X(X)  = 


m 


»(X) 


for  m=l:M  , 


3X„ 

m  m 


X(X) 


//__( X) 
mm 


//(X) 


//m(M  VX) 


//(X)  n(\) 


for  m,m=l:M 


(B-3) 


The  MxM  HM  of  joint  CGF  x(X)  is,  for  X  s  ROA(//), 


a^  1 

ax  ax 

m  m 

//(X)  ■ 

//(X)  tj{\) 

B-1 


There  follows,  from  equations  (B-l)  and  {B-2), 

~  J  ^m  z)  p{2)  for  m=l:M  , 


*  J  ^m  z)  p(z)  for  m,ni=l:M 


(B-5) 


All  these  integrals  converge  because  X  e  R 

r  u' 


Define  (tilted  PDF)  function 


p(z,X)  = 


=  ®xp(X^  2)  Dfz) 


U{\) 


for  X  e  ROA{//)  , 


and  define  the  quantities 


(B-6) 


“  J  ^m  m=l:M  . 


(B-7) 


Then,  the  m,m  element  of  matrix  H(X)  is,  from  equation  (B-3), 


“5i(X)  “  J  ^m  ^m  ~  m,m=l;M 


(B-8) 


and  X  e  ROA(//);  that  is,  X^.  =  Re{X)  e  R^. 

Now,  let  vector  X  be  real;  that  is,  X.  -  0  and  X  =  X^.  e  r  . 
It  is  presumed  that  there  are  no  linear  dependencies  among 
components  {2^}  of  RV  2.  Also,  let  real  Mxl  vectors 


A/(Xj,)  =  t^^(Xj,)  ...  ;/„(Xj.)]^  ,  a  =  [a^ 


V  =  2  -  ri{\)  . 


Then,  linear  combination 


(B-9) 


B-2 


f 


with  probability  1  . 


(B-10) 


T 

a  V  7^  0  for  any  a  ^  0 
Then,  the  following  average  must  be  positive; 


A(Xj.) 


/  T  x2 

( a  v) 


exp[x^  z] 

//( Xj.) 


>  0 


for  any  a  7^  0 


(B-ll) 


because  the  exp  term  and  are  always  positive.  Upon 

expansion, 

A{Xj.)  =  a"^  e[v  v"^  exp^X^  zj/t/CXj.)]  a  s  a^  C(Xj.)  a  >  0  (B-12) 

for  any  a  7^  0.  Therefore,  MxM  matrix  C(X  )  is  PD  for  X  c  R  . 

r  z  fj 


The  m,m  element  of  C(Xj,)  is,  using  equations  (B-9)  and  (B-6), 

■  J  K  -  ')  P(^)A<X,)  - 

.  J  dz  [z^  -  J„(X^)][zj,  -  2„(Xj.)]  g(z,X^)  . 

■  J  ■*'  K  'm  -  2m'\>  =5  -  %<\»]  * 

“  J  -  2„(Xj)  ?ji(Xj.)  for  .  (B-13) 

It  follows  from  equation  (B-8)  that 

m,m=l:M  ,  (B-14) 

meaning  that  H(Xj^)  *  C(Xj,)  is  a  PD  matrix.  That  is,  HM  H(X)  of 
joint  CGF  x(X)  is  PD  for  X  real  and  in  R  . 
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APPENDIX  C  —  POSITIVE  CURVATURE  OP  LOGARITHM  OF  INTEGRAND 
Consider  the  one-dimensional  Gaussian  RV  with  MGF  and  CGF 

/j{\)  =  exp^jff^X^  +  mxj  , 

X(X)  =  In  A/(X)  =  Iff^X^  +  mX  ,  (Cwi) 

for  all  X;  that  is,  Then,  for  evaluation  of  the  PDF 

at  field  point  2,  the  logarithm  A(X,2)  of  the  integrand  Y(X,z) 
and  its  derivatives  (with  respect  to  X)  are 

A(Xj.,z)  =  X{\)  -  zXj.  =  +  (m-z)Xj,  , 

A'(X^,z)  =  X'(Xj,)  -  z  =  ™  “  z  , 

2 

A”(X^^z)  ==  X"(X^)  =  a  ,  independent  of  z  ,  (C-2) 

for  all  X^ •  Observe  that  A"(X^/Z)  >  0  for  all  X^,  provided  that 
2 

a  >  0;  that  is,  the  curvature  of  A(X^,z)  is  positive  for  all  X^. 
Therefore,  A(X^,z)  can  only  have  a  single  minimum  on  the  real  X 
axis,  which  occurs  at  X(z)  =  (z-m)/<y^. 

For  the  next  example,  consider  the  one-dimensional  PDF 
exp(— z)  for  z  >  0*  The  pertinent  X-domain  functions  are 

/^(X^)  =  1/(1^X^)  ,  A(X^,z)  =  -ln(l-X^)  -  zX^  , 

A'(Xj.,2)  =  l/d-Xj.)  -  z  ,  A"(Xj.,z)  -  l/d-Xj.)^  ,  (C-3) 

for  Xj.  <  1;  that  is,  =  (-",1).  Again,  curvature  A"{Xj.,z)  >  0 


C-1 


for  all  Xj.  <  1,  although  the  curvature  tends  to  0+  as  ^  -®. 
Therefore,  A(Xj^,z)  can,  at  most,  have  only  one  minimum  for 

<  1,  but  it  may  have  none.  For  example,  if  z  >  0,  the  minimum 
of  A(Xj.,z)  is  at  X(z)  =  1  -  l/z;  as  z  0+,  then  X(z)  -®. 

On  the  other  hand,  if  z  <  0,  A{\^,z)  is  a  strictly 

monotonically  increasing  function  in  R^,  such  that  A(-®,z)  =  -®, 
A(l-,z)  =  +•.  This  last  region  for  z  (<  0)  corresponds  to  field 
points  where  the  exponential  PDF  is  zero;  that  is,  z  <  0  is  an 

unfeasible  point  for  the  PDF,  or  a  point  where  the  PDF  is  zero. 


Both  examples  above  illustrate  that  if  A’'(Xj,,z)  >  0  for  all 
Xr  e  R^,  then  A(Xj.,z)  can  have  only  one  minimum,  at  most,  in  that 

region;  that  is,  in  the  neighborhood  of  any  point  X  c  R 

a  fj' 

M\,z)  -  A’(X^,Z)  (X^-X^)  +  fA"(X^,2;)  (X_.-X,)^  .  (c-4) 

The  positiveness  of  the  last  coefficient  for  all  X^^  e  R 

indicates  that  the  function  A(X^,z)  has  positive  curvature  at  all 
points  in  R^. 


in  M  dimensions,  with  X  -  [ X^  ...  X^]^,  equation  (c-4) 
generalizes  to 


A(Xj.,z)  £  A(X  ,z)  +  G(X^,z)'^(X  -X  )+4 


H(Xg,z)(Xj.-X^)  (C-5) 


for  X^  e  R^,  where  Mxl  gradient  vector 


G(X,z ) 


(C-6) 


and  MxM  HM 


H(X,z) 


HI  m 


A( X,z) 


m ,  iii=l :  M  . 


(C-7) 


If  matrix  H(X^,z)  is  PD,  the  quadratic  form  in  equation  (C-5)  is 
positive  for  all  jt  X^;  that  is,  A(X^,z)  has  positive  curvature 
at  X  =  X  .  If  matrix  H(X  ,z)  is  PD  for  all  X  c  R  ,  then 

^  Cl  L  jT  fj 

A(X^,z)  has  positive  curvature  at  all  X^  e  R^,  and  can  have,  at 
most,  one  minimum  in  this  MD  region  R^.  Again,  this  minimum  can 
occur  for  some  components  of  X^,  equal  to  infinity,  or  some 
components  of  gradient  G(Xj,,z)  may  be  nonzero.  These  latter  two 
cases  correspond,  respectively,  to  the  situation  where  the 
statistical  average  a  of  interest  is  zero,  or  the  integrand 
Y(Xj,,z)  has  no  real  SP  for  X^,  e  R^. 
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APPENDIX  D  -  COMPLEX  LOCATIONS  OF  SADDLEPOINTS 


Consider  PDF  p(z)  =  exp(-z)  [1  +  cos(az)]  for  z  >  0.  This 
PDF  is  nonnegative  for  z  >  0;  although  its  area  is  not  unity, 
scaling  is  irrelevant  here.  Let  y  =  1  -  X  and  b  =  a^.  Then  the 
corresponding  MGF  and  CGF  and  their  derivatives  are,  for  <  1, 

,,X)  .  ^4^  ,  ^.(X)  -  .  «"(X)  -  , 

y(y^+b)  y^(y%b)2  y^y^+b)^ 


XMX) 


//(X)  ■ 


2v^+bv^+b^ 

2  2  ' 


y(y‘‘+b)(2y‘'+b) 


X"  (  X) 


4v^+7v^b^+8v^b^+b^ 

y^(y^+b)^(2y2+b)^ 


(D-1) 


Consider  the  ACs  of  all  these  functions,  and  observe  that 
X"(Xj,)  >  0  for  all  X^. 


The  corresponding  properties  for  the  integrand  and  its 
logarithm  are  given  by 

Y(X,z)  =  )i/(X)  exp(-Xz)  , 

A(X,z)  =  In  Y(X,z)  -  x(X)  -  Xz  , 

A'(X,z)  =  X'(X)  -  z  ,  A"(X)  =  X"(X)  .  (D-2) 

In  order  to  find  the  SPs  of  integrand  Y(X,z),  set  A'(X,2)  =  0  or 

A 

X'(X)  =  2  >  0,  and  obtain  the  equation 

2zy^  -  2y^  +  3zby^  -  by^  +  zb^y  -  b^  =  0  ,  (D-3) 


D-1 


with  X  =  1  -  y. 

As  an  example,  consider  a  =  1  and  z  =  0.1.  Equation  (D-3) 
then  yields  the  five  SPs 

X  =  -8.8995  ,  0.4951  ±  i  0.70305  ,  1.4547  ±  i  0.68369  .  (D-4) 

There  is  one  real  root  and  one  complex  pair  of  roots  within  the 
R0A(//);  that  IS,  with  X^.  <  1.  There  is  another  complex  pair  with 
Xj.  >  1,  which  is  outside  the  original  ROA(//).  The  SP  at  -8.8995 
connects  -i«  to  +i»  directly,  by  means  of  a  BC.  Use  of  the  other 
four  SPs  requires  using  the  AC  of  yt/(X)  and  the  valley  at  X  =  +-, 
in  addition  to  the  residue  of  the  pole  of  Y(X,z)  at  X  =  1. 
Alternatively,  the  contour  could  be  moved  to  X  =  +»  and  the  three 
residues  added;  the  result  would  be  the  exact  PDF  p(z),  with  no 
need  to  resort  to  SPs  or  a  SPA  at  all. 


An  example  of  a  PDF  with  additional  real  SPs  outside  the 
ROA(^)  is  furnished  by  the  sum  of  three  weighted,  independent 
unit-variance  exponential  RVs;  namely,  z  -  e^^  +  02/2  +  63/3. 
The  PDF  and  MGF  of  RV  2  are 

P(z)  -  3  exp(-z)  (1  -  exp(-z)]2  for  z  >  0  , 

~  (l-X)(l-X/2)(l-X/3)  Xj,  <  1  . 

There  follows 


AMX,z) 


z 


+ 


1 

1-X 


+ 


1 

2:^ 


1 

3-X 


(D-6) 


D-2 


The  numerator  of  A'(X,z)  is 
3  2 

X  z  +  X  (3  -  6  z)  +  X  (11  z  -  12)  +  (11  -  6  z)  .  (D-7) 

For  example,  with  z  =  13/12,  the  numerator  in  equation  (D-7)  has 
zeros  at  X  =  -1  and  at  X  =  (55  ±  /217)/26  *  1.54881  and  2.68196, 
all  of  which  are  on  the  real  X-axis.  The  latter  two  SPs  are 
outside  the  ROA(/j),  which  is  X^.  <  1.  These  SPs  are  the  only  ones 
for  integrand  T(X,z);  there  are  no  complex  SPs  in  this  example. 

If  desired,  the  original  BC  in  the  ROA(/u)  could  be  moved  so 
as  to  pass  through  the  SP  of  Y(X,z)  at  X  =  1.54881  when  the 
residue  of  the  AC  of  Y(X,z)  at  pole  location  X  =  1  is  also 
accounted  for.  Alternatively,  the  BC  could  be  moved  further  to 
the  right  so  as  to  pass  through  the  SP  of  Y(X,z)  at  X  =  2.68196 
when  the  residues  at  both  X  =  1  and  X  -  2  are  also  accounted  for. 
Finally,  the  BC  could  be  moved  off  to  the  valley  of  Y(X,z)  at 
X  =  +®,  giving  an  exact  result  for  the  PDF  in  terms  of  the  three 
residues  of  the  AC  of  Y(X,z)  at  pole  locations  X  *  1,  2,  and  3. 
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APPENDIX  E  -  ABSENCE  OF  SADDLEPOINT 
EVEN  THOUGH  z  IS  A  FEASIBLE  FIELD  POINT 


Consider  the  one-dimensional  PDF  and  MGF  pair: 


exp(-z ) 
(1+z)” 


for  z  >  0  ; 


(E-1) 


"  J  “  exp(w)  Ej^(w)  ,  w  =  1-X  (  e-2  ) 

upon  use  of  reference  3,  equation  5.1.4.  The  integral  in 
equation  (E-2)  converges  for  X^.  <  1  for  m  >  1,  and  for  X^.  <  1  for 
m  <  1.  (Strictly  speaking,  does  not  have  unit  area;  its 

area  is  /^jj^(O)  =  e  Ej^(l),  However,  absolute  scaling  is  not 
important  to  the  concepts  of  this  appendix.) 


A  useful  recursion  is 

^m^^^  “  [l  -  w  />i„_i(X)]  for  m  >  1  ,  (E-3) 

along  with  starting  values 

^0^^^  “  w  '  ^i^X)  =  exp(ft))  E^(w)  ,  both  for  «^  >  0  .  (E-4) 

In  particular,  for  >  0, 

//2(X)  =  1  -  w  exp(a))  Ej^(w)  ,  (E-5) 

^3^^^  “  -  ft)  +  «  exp(ft))  E^(ft))j  .  (E-6) 
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The  ACs  of  all  these  functions  behave  as  1/(1-X)  as  X 
according  to  reference  3,  equation  5.1.51. 


Additional  useful  relations  are 

' 

for  m  >  1  , 

*  Ur-aHm-l)  £or  m  >  2  . 


(E-7) 


{E-8) 


(E-9) 


The  derivative  of  the  corresponding  CGF  is 
V-,1,  1  , 

^  m  >  2  .  (E-10) 


Therefore,  z  =  l/(m-2)  is  a  breakpoint  as  to  whether  a  real  SP  of 
integrand  //^(X)  exp(-zX)  can  exist  for  m  >  2.  In  addition. 


m-l 

(m-3 ) (m-2)^ 


for  m  >  3  , 


(E-11) 


_ 2 _ 

(m-3 ) (m-2 ) (m-l ) 


for  m  >  3  . 


(E-12) 


For  m  .  2,  any  z  >  0  Is  allowed:  a  real  SP  always  exists. 
That  Is,  integrand  fiX^.z)  =  exp(-sX_.)  always  has  a 

minimum  for  X^.  <  1,  when  z  >  0. 

From  now  on,  let  m  >  2.  Then,  there  can  be  no  real  SP  of 


E-2 


integrand  T(X,z)  —  //jjj(X)  exp(— zX)  if  z  >  l/(m— 2)  because  of 
relation  (E-10).  That  is,  the  curve  of  Y(Xj,,z)  versus  X^  is 
still  monotonically  decreasing  at  X^  =  1  if  z  >  l/(m-2).  This 
feature  is  illustrated  in  figure  E-1  for  m  =  3.  The  slope  of 
Y(Xj.,z)  at  Xj.  =  1  is,  by  use  of  equations  (E-8)  and  {E-9), 

exp(-z)  r  1  1  ^ 

m-1  Li^  "  m  >  2  .  (E-13) 

Then,  the  minimum  of  the  integrand  (but  not  a  zero  slope)  is 
realized  at  the  boundary  of  definition,  namely,  at  =  1.  The 
point  X  =  1  is  a  BP  of  the  AC  of  yWjjj(X)  because  Ej^(w)  has  a 
logarithmic  singularity  at  w  =  0;  that  is,  at  X  =  1.  The  product 
(0  E^(oo)  =  0  at  «  =  0. 

An  approach  similar  to  appendix  A  can  be  employed  here.  As 
an  example,  consider  m  =  3,  as  in  equation  (E-6).  The  AC  of 
//^(X)  has  a  logarithmic  BP  at  X  *  1  according  to  reference  3, 
equation  5.1.11.  Let  z  >  0  and  move  the  BC  to  a  keyhole  contour 
wrapped  around  X  =  1  and  the  positive  real  axis;  this  is  the 
steepest  descent  contour.  Then,  since  ~  1/(1-X)  as  X  ®, 

the  integrals  on  the  large  circular  arcs  tend  to  zero  as  the 
radius  tends  to  infinity.  Although  Ej^(w)  has  a  logarithmic 
singularity  at  X  =  1,  the  quantity  (1-X)^  tends  to  zero  faster, 
so  that  the  integral  on  the  small  circle  around  X  =  1  tends  to 
zero  as  the  radius  tends  to  zero.  On  the  upper  (lower) 
horizontal  line,  X^.  >  i,  E^(l-X)  =  -  Ei(X-l)  +(-)  in.  Since  the 
integrals  are  in  opposite  directions,  only  the  imaginary  parts 
remain,  giving  the  exact  result 


E-3 


^3^^^  “  2  J  dX  exp(-2X)  exp(l-X)  in  »  e2E£iz:|i  . 


(1  +  2) 


The  alternative  two-sided  PDF  £^(2)  has  a  somewhat  differs 
behavior : 


D  (2)  = 

(l+|z|)® 


for  all  2  . 


(E-15: 


The  MGP  iu^(\)  exists  only  for  -1  <  X^  <  1.  m  fact, 

-1  <  <  1  . 


(E-16) 


For  m  -  3,  as  X^.  varies  in  t-1,1],  integrand  £jn(Xj.)  exp(-2X^) 

varies  over  a  f^tnite  range;  it  does  not  realise  a  minimum  when 

|z|  >  0.61.  Therefore,  no  real  SP  exists  for  this  integrand  when 
jz I  >  0.61. 


'F{Xr,Z)  0.15 


zero  slope  ■ 


0  0.1  0.2  0.3  0.4  0.5  0.6  0.7  ; 

Xr 

Figure  E-1.  Integrand  of  Inverse  Laplace  Transform  for  m  . 


APPENDIX  F  -  HESSIAN  MATRIX  OF  In  y(X)  IS 
NONNEGATIVE  DEFINITE  FOR  X  REAL  AND  IN  R 

Y 

Let  g(u)  be  a  real  nonnegative  scalar  function  of  MD  vector 
u;  g(u)  need  not  have  finite  volume.  Define  MD  X  =  X^.  +  i  and 

y(X)  »  J  du  exp(-x'^  u)  g(u)  .  (F_l) 


Let  this  MD  integral  converge  for  X^,  e  R^.  Define 


y(X)  =  Y_(X)  for  m=l;M  . 
m 


(F-2) 


There  follows 


^  In  y(X) 
ni 


y(X) 


for  m=l;M  , 


(F-3) 


and 


ax^  8X 
HI  in 


y(X) 


TTxT  TIxT  :M  .  (P-4) 


The  MxM  HM  of  In  y(X)  is,  for  X  e  R  , 

r  y' 


H(X)  H 


3X„  ax 

in  m 


y(X)  y(X)  y(X) 


(F-5) 


From  equations  (F-1)  and  (F-2),  there  follows 

^m^^^  “  “  J  %  exp(-x'*^  u)  g(u)  for  m=l:M  (F-6) 

and 

’’^mm^^^  “  J  ^m  %  exp(-X^  u)  g(u)  for  m,m=l:M  .  (F-7) 


F-1 


Now,  for  Re(X)  e  R  ,  define  (tilted)  functi( 


g(u,X)  =  exp(-X  u)  Qfu) 
y(  X) 

and  the  quantities 


(F-8) 


“  J  %  9(u,X)  for  m=l:M  . 

Then,  the  in,m  element  of  HM  H{X)  in  equation  (F-5)  becomes 
®mm^^^  ~  J  ’^m  %  9(u,X)  -  v  (X)  v  (X)  for  m,m=l:M  . 


(F-9) 


(F-10) 


Observe  that  the  function  g(u,X)  defined  in  equation  (f-8) 
has  unit  volume  in  u  for  X  =  X^  e  R  ;  that  is. 


J  du  g(u,X^)  =  1  for  all  X^ 


e  R  . 
Y 


(F-11) 


Thus,  §(u,X^)  has  the  properties  of  a  legal  PDF  of  arguaient  u  for 
any  s  E^,  in  particular,  g(u,X^)  is  real,  nonnegative,  and  has 
unit  volume  in  «d  u  space.  Therefore,  considering  §(u,X^)  as  a 
MD  PDF  in  u,  the  mean  of  the  corresponding  m-th  Rv  u^(X  )  is 


|(Xr)  -  J  du  u^  g{u,X^)  for  m»l:n  , 


(F-12) 


while  the  covariance  of  the  m,m  pair  of 


RVs  u^(X^),  u^(x^)  is 


~mm^  r^  J  m,m=l:M. 

But,  the  covariance  matrix  of  any  pdf  is  always  nonnegative 
definite  because  with  random  vector  w  =  u(  X^. )  -  v(  X^. )  = 


(F-13) 


F-2 


(F-14) 


***  average 

E  (a*^  w)^  =  a'^  E(w  w*^)  a  =  a^  U(Xj.)  a  >  0 

for  any  Mxl  real  vector  a.  U(Xj.)  = 
covariance  matrix  of  RV  u(X^)  and  must  be  nonnegative  definite. 

Comparison  of  equations  (F-10)  and  (F-13)  reveals  that 

m,m-l:M  .  (F-15) 

Therefore,  matrix  H(Xj,)  =  U(X^)  is  nonnegative  definite;  that  is, 

the  HM  H(X)  of  In  r{\)  is  nonnegative  definite  for  X  real  and  in 

R  . 

r 


As  a  first  example,  consider 

g(u)  =  5{u^  -  2^)  ...  8(Uj^  -  z^)  ,  (F-16) 

for  which 

y(X)  =  exp(-X^  X)  for  all  X  .  (F-17) 

The  HM  of  In  y(X)  (=  -X*^  2)  is  2ero  for  all  X;  this  HM  has  M  2ero 
eigenvalues  and  is  nonnegative  definite  for  all  X. 

The  second  example  is 

g(u)  =  5(Uj^  -  z^)  U(U2  -  Z2)  •..  U(Ujj  -  z^)  ,  (F-18) 

for  which 


F-3 


(F-19) 


r(X)  =  exp(-X^  >  o  ,  m=2:M  , 


and  X^  is  arbitrary.  There  follows 


M 


In  y(X)  =  -X  2  -  II]  , 


in=2 


for  which  the  gradient  vector  ii 


G(X)  -  [-2^  -22-VX2  •••  -2„-l/Xj^]^ 


The  MxM  HM  of  In  Y{X)  is  then 


H(X)  =  diag[o  1/X^  ...  i/x2j  , 

The  matrix  H(X^)  has  one  zero  eigenvalue  and  M-1  positive 
eigenvalues;  thus,  H(Xj.)  is  nonnegative  definite. 


(F-20) 


{F-21) 


F-22) 
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